/* Copyright 2012 Tobias Marschall
 *
 * This file is part of CLEVER.
 *
 * CLEVER is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CLEVER is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CLEVER.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <cassert>
#include <algorithm>

#include "SplitAligner.h"

using namespace std;

SplitAligner::SplitAligner(int cost_threshold, gap_costs_t gap_costs, int indel_cost, int gap_extend_cost, int split_cost, int interchrom_split_cost, int min_anchor_length, int secondary_aln_phred_diff, int phred_offset, bool simple_cigars, bool rightmost, bool allow_soft_clipping, int soft_clip_open_cost, int soft_clip_extend_cost) {
	assert(secondary_aln_phred_diff >= 0);
	assert((int)std::numeric_limits<cost_t>::max() > cost_threshold);
	this->cost_threshold = cost_threshold;
	this->gap_costs = gap_costs;
	this->indel_cost = indel_cost;
	this->gap_extend_cost = gap_extend_cost;
	this->split_cost = split_cost;
	this->interchrom_split_cost = interchrom_split_cost;
	this->min_anchor_length = min_anchor_length;
	this->phred_offset = phred_offset;
	this->secondary_aln_phred_diff = secondary_aln_phred_diff;
	this->simple_cigars = simple_cigars;
	this->rightmost = rightmost;
	this->allow_soft_clipping = allow_soft_clipping;
	this->soft_clip_open_cost = soft_clip_open_cost;
	this->soft_clip_extend_cost = soft_clip_extend_cost;
	switch (gap_costs) {
		case LINEAR:
			this->max_indels = cost_threshold / indel_cost;
			break;
		case AFFINE:
			this->max_indels = max(0, (cost_threshold - indel_cost) / gap_extend_cost);
			break;
		default:
			assert(false);
	}
}

SplitAligner::~SplitAligner() {
}

auto_ptr<vector<SplitAligner::split_alignment_t> > SplitAligner::combineAnchors(const anchor_alignment_t& left_anchor, const anchor_alignment_t& right_anchor) const {
	assert(left_anchor.total_query_length == right_anchor.total_query_length);
	int length = left_anchor.total_query_length;
	auto_ptr<vector<SplitAligner::split_alignment_t> > result(new vector<SplitAligner::split_alignment_t>());
	// =============== PART 1: determine deletion split reads ============================
	// determine start and end of region of overlap of columns
	int start_column = max((int)(length - right_anchor.column_count + 1), min_anchor_length);
	int end_column = min((int)(left_anchor.column_count - 1), length - min_anchor_length);
	assert(start_column >= 0);
	assert(end_column <= length);
	if (start_column <= end_column) {
		int right_col = length - start_column;
		int best_score = -1;
		for (int left_col = start_column; left_col <= end_column; ++left_col, --right_col) {
			if (best_score!=-1) {
				if (left_anchor.column_minima[left_col] + right_anchor.column_minima[right_col] + split_cost > min(best_score + secondary_aln_phred_diff, cost_threshold) ) continue;
			} else {
				if (left_anchor.column_minima[left_col] + right_anchor.column_minima[right_col] + split_cost > cost_threshold) continue;
			}
			// cerr << "left_col: " << left_col << ", right_col: " << right_col << endl;
			for (int left_rel_row = -max_indels; left_rel_row <= max_indels; ++left_rel_row) {
				// the diagonal in the dp matrix
				int left_diag = left_rel_row + max_indels;
				int left_cost = left_anchor.diagonals[left_diag][left_col];
				if (left_cost > left_anchor.column_minima[left_col] + secondary_aln_phred_diff) continue;
				int left_diff = left_cost - left_anchor.column_minima[left_col];
				// the current reference position on the left anchor (incl.)
				int left_ref_pos = left_anchor.ref_start_pos + left_col + left_rel_row - 1;
				for (int right_rel_row = -max_indels; right_rel_row <= max_indels; ++right_rel_row) {
					int right_diag = right_rel_row + max_indels;
					int right_cost = right_anchor.diagonals[right_diag][right_col];
					if (right_cost > right_anchor.column_minima[right_col] + secondary_aln_phred_diff - left_diff) continue;
					// the current reference position on the right anchor (incl.)
					int right_ref_pos = right_anchor.ref_start_pos - (right_col + right_rel_row - 1);
					// cerr << "  left_ref_pos: " << left_ref_pos << ", right_ref_pos: " << right_ref_pos << endl;
					if (left_ref_pos + 1 >= right_ref_pos) continue;
					// total cost for the split alignment
					int total_cost = left_cost + right_cost + split_cost;
					if (best_score != -1) {
						if (total_cost <= best_score + secondary_aln_phred_diff) {
							push_deletion(left_anchor, right_anchor,  left_ref_pos, right_ref_pos, left_diag, right_diag, left_col, right_col, total_cost, result.get());
						}
					} else {
						if (total_cost < cost_threshold) {
							push_deletion(left_anchor, right_anchor,  left_ref_pos, right_ref_pos, left_diag, right_diag, left_col, right_col, total_cost, result.get());
						}
					}
				}
			}
		}
	}
	// =============== PART 2: determine insertion split reads ============================
	// determine region of possible breakpoint positions in the reference
	int start_ref_pos = max(right_anchor.ref_start_pos - right_anchor.row_count, left_anchor.ref_start_pos);
	int end_ref_pos = min(left_anchor.ref_start_pos + left_anchor.row_count - 1, right_anchor.ref_start_pos - 1);
	// cerr << "start_ref_pos: " << start_ref_pos << ", end_ref_pos: " << end_ref_pos << endl;
	if (start_ref_pos <= end_ref_pos) {
		int ref_pos = start_ref_pos;
		int left_row_idx = ref_pos - left_anchor.ref_start_pos + 1;
		int right_row_idx = right_anchor.ref_start_pos - ref_pos;
		int best_score = -1;
		for (; ref_pos <= end_ref_pos; ++ref_pos, ++left_row_idx, --right_row_idx) {
			if (best_score != -1) {
				if (left_anchor.row_minima[left_row_idx] + right_anchor.row_minima[right_row_idx] + split_cost > min(best_score + secondary_aln_phred_diff, cost_threshold) ) continue;
			} else {
				if (left_anchor.row_minima[left_row_idx] + right_anchor.row_minima[right_row_idx] + split_cost > cost_threshold) continue;
			}
			// cerr << "left_row_idx: " << left_row_idx << ", right_row_idx: " << right_row_idx << endl;
			for (int left_diag_idx = 0; left_diag_idx < left_anchor.diagonals.size(); ++left_diag_idx) {
				const std::vector<cost_t>& left_d = left_anchor.diagonals[left_diag_idx];
				int left_col_idx = left_row_idx - left_diag_idx + max_indels;
				if ((left_col_idx < 0) || (left_col_idx >= left_d.size())) continue;
				if (left_col_idx < min_anchor_length) continue;
				int left_cost = left_d[left_col_idx];
				if (left_cost > left_anchor.row_minima[left_row_idx] + secondary_aln_phred_diff) continue;
				int left_diff = left_cost - left_anchor.row_minima[left_row_idx];
				int left_query_pos = left_col_idx - 1;
				for (int right_diag_idx = 0; right_diag_idx < right_anchor.diagonals.size(); ++right_diag_idx) {
					const std::vector<cost_t>& right_d = right_anchor.diagonals[right_diag_idx];
					// diag = row - column + max_indels;
					int right_col_idx = right_row_idx - right_diag_idx + max_indels;
					if (right_col_idx < min_anchor_length) continue;
					if ((right_col_idx < 0) || (right_col_idx >= right_d.size())) continue;
					int right_query_pos = right_anchor.total_query_length - right_col_idx;
					if (left_query_pos + 1 >= right_query_pos) continue;
					int right_cost = right_d[right_col_idx];
					if (right_cost > right_anchor.row_minima[right_row_idx] + secondary_aln_phred_diff - left_diff) continue;
					int total_cost = left_cost + right_cost + split_cost;
					if (best_score != -1) {
						if (total_cost <= best_score + secondary_aln_phred_diff) {
							push_insertion(left_anchor, right_anchor, ref_pos + 1, left_query_pos, right_query_pos, left_diag_idx, right_diag_idx, left_col_idx, right_col_idx, total_cost, result.get());
						}
					} else {
						if (total_cost < cost_threshold) {
							push_insertion(left_anchor, right_anchor, ref_pos + 1, left_query_pos, right_query_pos, left_diag_idx, right_diag_idx, left_col_idx, right_col_idx, total_cost, result.get());
						}
					}
				}
			}
		}
	}
	// =============== PART 3: sort results and only retain sufficiently good alignments ============================
	if (result->size() > 0) {
		sort(result->begin(), result->end(), split_alignment_comparator_t(rightmost));
		if (secondary_aln_phred_diff != -1) {
			vector<SplitAligner::split_alignment_t>::iterator it = result->begin();
			for (; it != result->end(); ++it) {
				if (it->phred_score > result->at(0).phred_score + secondary_aln_phred_diff) {
					result->erase(it, result->end());
					break;
				}
			}
		}
	}
	return result;
}

bool SplitAligner::doInterchromosomalSearch(int best_known_score) const {
	return interchrom_split_cost - secondary_aln_phred_diff <= best_known_score;
}

auto_ptr<vector<SplitAligner::interchrom_split_alignment_t> > SplitAligner::combineAnchorsInterchromosomal(const anchor_alignment_t& anchor1, const anchor_alignment_t& anchor2, int best_known_score) const {
// 	cerr << "SplitAligner::combineAnchorsInterchromosomal, " << interchrom_split_cost << endl;
	assert(anchor1.total_query_length == anchor2.total_query_length);
	int length = anchor1.total_query_length;
	auto_ptr<vector<SplitAligner::interchrom_split_alignment_t> > result(new vector<SplitAligner::interchrom_split_alignment_t>());
	// determine start and end of region of overlap of columns
	int start_column = max((int)(length - anchor2.column_count + 1), min_anchor_length);
	int end_column = min((int)(anchor1.column_count - 1), length - min_anchor_length);
	assert(start_column >= 0);
	assert(end_column <= length);
	if (start_column <= end_column) {
		int col2 = length - start_column;
		int best_score = best_known_score;
// 		cerr << "best_known_score: " << best_known_score << endl;
		for (int col1 = start_column; col1 <= end_column; ++col1, --col2) {
			if ((best_score!=-1) && (secondary_aln_phred_diff!=-1)) {
				if (anchor1.column_minima[col1] + anchor2.column_minima[col2] + interchrom_split_cost > min(best_score + secondary_aln_phred_diff, cost_threshold) ) continue;
			} else {
				if (anchor1.column_minima[col1] + anchor2.column_minima[col2] + interchrom_split_cost > cost_threshold) continue;
			}
// 			cerr << "left_col: " << col1 << ", right_col: " << col2 << endl;
			for (int rel_row1 = -max_indels; rel_row1 <= max_indels; ++rel_row1) {
				// the current reference position on the left anchor (incl.)
				int breakpoint1_ref_pos = anchor1.ref_start_pos + (anchor1.type==LEFT_ANCHOR?1:-1) * (col1 + rel_row1 - 1);
				// the diagonal in the dp matrix
				int diag1 = rel_row1 + max_indels;
				int cost1 = anchor1.diagonals[diag1][col1];
				for (int rel_row2 = -max_indels; rel_row2 <= max_indels; ++rel_row2) {
					// TODO: Speedup: no need to examine all combinations
					// the current reference position on the right anchor (incl.)
					int breakpoint2_ref_pos = anchor2.ref_start_pos + (anchor2.type==LEFT_ANCHOR?1:-1) * (col2 + rel_row2 - 1);
					int diag2 = rel_row2 + max_indels;
					int cost2 = anchor2.diagonals[diag2][col2];
					// if (left_ref_pos + 1 >= right_ref_pos) continue;
					// total cost for the split alignment
					int total_cost = cost1 + cost2 + interchrom_split_cost;
// 					cerr << "  left_ref_pos: " << breakpoint1_ref_pos << ", right_ref_pos: " << breakpoint2_ref_pos << " --> " << total_cost << endl;
					if ((best_score!=-1) && (secondary_aln_phred_diff!=-1)) {
						if (total_cost <= min(best_score + secondary_aln_phred_diff, cost_threshold)) {
							push_interchromosomal(anchor1, anchor2,  breakpoint1_ref_pos, breakpoint2_ref_pos, diag1, diag2, col1, col2, total_cost, result.get());
						}
					} else {
						if (total_cost < cost_threshold) {
							push_interchromosomal(anchor1, anchor2,  breakpoint1_ref_pos, breakpoint2_ref_pos, diag1, diag2, col1, col2, total_cost, result.get());
						}
					}
				}
			}
		}
	}
	// =============== sort results and only retain sufficiently good alignments ============================
	if (result->size() > 0) {
		sort(result->begin(), result->end(), interchrom_split_alignment_comparator_t());
		if (secondary_aln_phred_diff != -1) {
			vector<SplitAligner::interchrom_split_alignment_t>::iterator it = result->begin();
			for (; it != result->end(); ++it) {
				if (it->phred_score > result->at(0).phred_score + secondary_aln_phred_diff) {
					result->erase(it, result->end());
					break;
				}
			}
		}
	}
	return result;
}

void SplitAligner::compute_reverse_cigar(const anchor_alignment_t& anchor, int diag, int col, std::vector<BamTools::CigarOp>* target, int* mismatch_phred_score) const {
	assert(mismatch_phred_score!=0);
	BamTools::CigarOp current_op('\0', 0);
	assert(diag >= 0);
	assert(diag < (int)anchor.backtrace.size());
	assert(col >= 0);
	assert(col < (int)anchor.backtrace[diag].size());
	while ((col!=0) || (diag!=max_indels)) {
		// cerr << "  col: " << col << ", row: " << (diag+col-max_indels) << ", diag: " << diag << ", op: " << anchor.backtrace[diag][col] << ", indel_len: " << anchor.indel_length_table[diag][col] << endl;
		char next_cigar_char = '\0';
		int next_cigar_length = 1;
		int c = 0;
		switch (anchor.backtrace[diag][col]) {
		case LEFT:
			next_cigar_char = 'I';
			if (gap_costs == AFFINE) {
				next_cigar_length = anchor.indel_length_table[diag][col];
				assert(next_cigar_length>0);
				col -= next_cigar_length;
				diag += next_cigar_length;
			} else {
				col -= 1;
				diag += 1;
			}
			break;
		case ABOVE:
			next_cigar_char = 'D';
			if (gap_costs == AFFINE) {
				next_cigar_length = anchor.indel_length_table[diag][col];
				assert(next_cigar_length>0);
				diag -= next_cigar_length;
			} else {
				diag -= 1;
			}
			break;
		case ABOVE_LEFT:
			c = anchor.diagonals[diag][col] - anchor.diagonals[diag][col-1];
			if (simple_cigars) {
				next_cigar_char = 'M';
			} else {
				next_cigar_char = (c==0?'=':'X');
			}
			*mismatch_phred_score += c;
			col -= 1;
			break;
		default:
			assert(false);
		}
		if (next_cigar_char == current_op.Type) {
			current_op.Length += next_cigar_length;
		} else {
			if (current_op.Length > 0) target->push_back(current_op);
			current_op.Type = next_cigar_char;
			current_op.Length = next_cigar_length;
		}
	}
	if (current_op.Length > 0) target->push_back(current_op);
}

void SplitAligner::push_deletion(const anchor_alignment_t& left, const anchor_alignment_t& right,  int left_ref_pos, int right_ref_pos, int left_diag, int right_diag, int left_col, int right_col, cost_t cost, std::vector<split_alignment_t>* target) const {
	int mismatch_phred_score = 0;
	vector<BamTools::CigarOp> left_cigar;
	compute_reverse_cigar(left, left_diag, left_col, &left_cigar, &mismatch_phred_score);
	vector<BamTools::CigarOp> right_cigar;
	compute_reverse_cigar(right, right_diag, right_col, &right_cigar, &mismatch_phred_score);
	vector<BamTools::CigarOp> middle_cigar;
	middle_cigar.push_back(BamTools::CigarOp('D', right_ref_pos - left_ref_pos - 1));
	target->push_back(split_alignment_t(SPLIT_DELETION, right_ref_pos - left_ref_pos - 1, right_ref_pos, left_col + 1, cost, mismatch_phred_score));
	split_alignment_t& result = target->at(target->size() - 1);
	append_cigar(left_cigar.rbegin(), left_cigar.rend(), &result.cigar);
	append_cigar(middle_cigar.begin(), middle_cigar.end(), &result.cigar);
	append_cigar(right_cigar.begin(), right_cigar.end(), &result.cigar);
	// cerr << "Deletion at " << left_ref_pos + 2 << " - " << right_ref_pos << ", cost: " << cost << endl;
}

void SplitAligner::push_insertion(const anchor_alignment_t& left_anchor, const anchor_alignment_t& right_anchor, int ref_breakpoint, int left_query_pos, int right_query_pos, int left_diag_idx, int right_diag_idx, int left_col_idx, int right_col_idx, cost_t cost, std::vector<split_alignment_t>* target) const {
	int length = right_query_pos - left_query_pos - 1;
	int mismatch_phred_score = 0;
	vector<BamTools::CigarOp> left_cigar;
	compute_reverse_cigar(left_anchor, left_diag_idx, left_col_idx, &left_cigar, &mismatch_phred_score);
	vector<BamTools::CigarOp> right_cigar;
	compute_reverse_cigar(right_anchor, right_diag_idx, right_col_idx, &right_cigar, &mismatch_phred_score);
	vector<BamTools::CigarOp> middle_cigar;
	middle_cigar.push_back(BamTools::CigarOp('I', length));
	target->push_back(split_alignment_t(SPLIT_INSERTION, length, ref_breakpoint, left_col_idx + 1, cost, mismatch_phred_score));
	split_alignment_t& result = target->at(target->size() - 1);
	append_cigar(left_cigar.rbegin(), left_cigar.rend(), &result.cigar);
	append_cigar(middle_cigar.begin(), middle_cigar.end(), &result.cigar);
	append_cigar(right_cigar.begin(), right_cigar.end(), &result.cigar);
	// cerr << "Insertion at " << ref_breakpoint << ", length " << length << ", cost: " << cost << endl;
}

void SplitAligner::anchor_to_cigar(const anchor_alignment_t& anchor, int diag, int col, std::vector<BamTools::CigarOp>* cigar, int* mismatch_phred_score) const {
	if (anchor.type == LEFT_ANCHOR) {
		vector<BamTools::CigarOp> rev_cigar;
		compute_reverse_cigar(anchor, diag, col, &rev_cigar, mismatch_phred_score);
		append_cigar(rev_cigar.rbegin(), rev_cigar.rend(), cigar);
	} else {
		compute_reverse_cigar(anchor, diag, col, cigar, mismatch_phred_score);
	}
}

void SplitAligner::push_interchromosomal(const anchor_alignment_t& anchor1, const anchor_alignment_t& anchor2,  int breakpoint1_ref_pos, int breakpoint2_ref_pos, int diag1, int diag2, int col1, int col2, cost_t cost, std::vector<interchrom_split_alignment_t>* target) const {
	assert(col1 + col2 == anchor1.total_query_length);
	assert(anchor1.total_query_length == anchor2.total_query_length);
	int mismatch_phred_score = 0;
	vector<BamTools::CigarOp> rev_cigar1;
	compute_reverse_cigar(anchor1, diag1, col1, &rev_cigar1, &mismatch_phred_score);
	vector<BamTools::CigarOp> soft_clip_cigar1;
	soft_clip_cigar1.push_back(BamTools::CigarOp('S', anchor1.total_query_length - col1));
	vector<BamTools::CigarOp> soft_clip_cigar2;
	soft_clip_cigar2.push_back(BamTools::CigarOp('S', anchor2.total_query_length - col2));
	vector<BamTools::CigarOp> rev_cigar2;
	compute_reverse_cigar(anchor2, diag2, col2, &rev_cigar2, &mismatch_phred_score);
	target->push_back(interchrom_split_alignment_t(breakpoint1_ref_pos, breakpoint2_ref_pos, cost, mismatch_phred_score, interchrom_split_cost));
	interchrom_split_alignment_t& result = target->at(target->size() - 1);
	if (anchor1.type == LEFT_ANCHOR) {
		append_cigar(rev_cigar1.rbegin(), rev_cigar1.rend(), &result.cigar1);
		append_cigar(soft_clip_cigar1.begin(), soft_clip_cigar1.end(), &result.cigar1);
	} else {
		append_cigar(soft_clip_cigar1.begin(), soft_clip_cigar1.end(), &result.cigar1);
		append_cigar(rev_cigar1.begin(), rev_cigar1.end(), &result.cigar1);
	}
	if (anchor2.type == LEFT_ANCHOR) {
		append_cigar(rev_cigar2.rbegin(), rev_cigar2.rend(), &result.cigar2);
		append_cigar(soft_clip_cigar2.begin(), soft_clip_cigar2.end(), &result.cigar2);
	} else {
		append_cigar(soft_clip_cigar2.begin(), soft_clip_cigar2.end(), &result.cigar2);
		append_cigar(rev_cigar2.begin(), rev_cigar2.end(), &result.cigar2);
	}
}

auto_ptr<SplitAligner::whole_alignment_t> SplitAligner::toWholeAlignment(const anchor_alignment_t& anchor) const {
	int length = anchor.total_query_length;
	// pick end columen
	int col = -1;
	bool soft_clipping = false;
	int soft_clipping_phred = 0;
	if (allow_soft_clipping) {
		// when soft clipping is allowed determine how much of the alignment to use
		int best_col = -1;
		int best_phred = this->cost_threshold + 1;
		for (int c = anchor.column_count - 1; c>=this->min_anchor_length; --c) {
			 int phred = anchor.column_minima[c];
			 if (c != length) {
				phred += soft_clip_open_cost + (length-c) * soft_clip_extend_cost;
			 }
			 if (phred < best_phred) {
				best_col = c;
				best_phred = phred;
			 }
		}
		if (best_col == -1) return auto_ptr<whole_alignment_t>(0);
		col = best_col;
		soft_clipping = col < length;
		if (soft_clipping) soft_clipping_phred = soft_clip_open_cost + (length-col) * soft_clip_extend_cost;
		// cerr << "Softclipping at " << col << " / " << length << " with cost " << best_phred << endl;
	} else {
		// if soft clipping is disabled, then use the full length
		if (anchor.column_count != length + 1) return auto_ptr<whole_alignment_t>(0);
		if (anchor.column_minima[length] > cost_threshold) return auto_ptr<whole_alignment_t>(0);
		col = length;
	}
	auto_ptr<whole_alignment_t> result(new whole_alignment_t());
	// pick best diagonal in chosen column of alignment. if two diagonals are equally good, pick
	// the one closest to the main diagonal
	int best_diag = -1;
	int best_phred = -1;
	for (int k=0; k<anchor.diagonals.size(); ++k) {
		int cost = anchor.diagonals[k][col];
		if (cost <= cost_threshold) {
			if ((best_diag==-1) || (cost < best_phred)) {
				best_diag = k;
				best_phred = cost;
			} else if ((cost==best_phred) && (abs(k-max_indels) < abs(best_diag-max_indels))) {
				best_diag = k;
				best_phred = cost;
			}
		}
	}
	// cerr << "  Selected diagonal " << best_diag << " with cost " << best_phred << " after adding soft_clip cost: " << (best_phred + soft_clipping_phred) <<  endl;
	if (soft_clipping) best_phred += soft_clipping_phred;
	assert(best_diag >= 0);
	assert((best_phred>=0) && (best_phred<=cost_threshold));
	result->phred_score = best_phred;
	vector<BamTools::CigarOp> cigar;
	int mismatch_phred_score = 0;
	compute_reverse_cigar(anchor, best_diag, col, &cigar, &mismatch_phred_score);
	result->mismatch_phred_score = mismatch_phred_score;
	if (anchor.type == LEFT_ANCHOR) {
		append_cigar(cigar.rbegin(), cigar.rend(), &result->cigar);
		if (soft_clipping) {
			BamTools::CigarOp clip_op('S', length - col);
			append_cigar(&clip_op, (&clip_op) + 1, &result->cigar);
		}
		result->start_position = anchor.ref_start_pos;
		result->end_position = anchor.ref_start_pos + col - 1 + (best_diag-max_indels);
	} else {
		assert(anchor.type == RIGHT_ANCHOR);
		if (soft_clipping) {
			BamTools::CigarOp clip_op('S', length - col);
			append_cigar(&clip_op, (&clip_op) + 1, &result->cigar);
		}
		append_cigar(cigar.begin(), cigar.end(), &result->cigar);
		result->start_position = anchor.ref_start_pos - col + 1 - (best_diag-max_indels);
		result->end_position = anchor.ref_start_pos;
	}
	return result;
}
